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ABSTRACT 



Numerical simulations suggest that the first galaxies are formed in protogalactic halos with virial temperatures > 10"^ K. It is likely that 
such halos are polluted with trace amounts of metals produced by the first generation of stars. The presence of dust can significantly 
change the chemistry and dynamics of early galaxies. In this article, we aim to assess the role of dust on the thermal and dynamical 
evolution of the first galaxies in the presence of a background UV flux, and its implications for the observability of Lyman alpha 
emitters and sub-mm sources. We have performed high resolution cosmological simulations using the adaptive mesh refinement code 
FLASH to accomplish this goal. We have developed a chemical network appropriate for these conditions and coupled it with the 
FLASH code. The main ingredients of our chemical model include the formation of molecules (both in the gas phase and on dust 
grains), a multi-level treatment of atomic hydrogen, line trapping of Lyman alpha photons and, photoionization and photodissociation 
processes in a UV background. We found that the formation of molecules (H2 and HD) is significantly enhanced in the presence of 
dust grains as compared to only gas phase reactions by up to two orders of magnitude. The presence of dust may thus establish a 
molecular ISM in high-redshift galaxies. The presence of a background UV flux strongly influences the formation of molecules by 
photodissociating them. We explore the evolution after a major merger, leading to the formation of a binary disk. These disks have gas 
masses of ~ 10^ M© at a redshift of 5.4. Each disk lies in a separate subhalo as a result of the merger event. The disks are supported 
by turbulent pressure due to the highly supersonic turbulence present in the halo. For values of J21 = 1000 (internal flux), we find 
that fragmentation may be enhanced due to thermal instabilities in the hot gas. The presence of dust does not significantly reduce the 
Lyman alpha emission. The emission of Lyman alpha is extended and originates from the envelope of the halo due to line trapping 
effects. We also find that dust masses of a few xlO^ M© are required to observe the dust continuum emission from z ~ 5 galaxies with 
ALMA. 
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1. Introduction 



Cosmology: theory - early Universe - Galaxies: formation - Atomic processes - Molecular 



The observation al frontier of high- r edshift galax ies is 

shifting rapidly (iBouwens et alJ 20101: iBunker et alJ l2010l: 
iBouwens et aDl201 lUVanzella et al.ll2'om . High-redshift galax- 
ies have been detected using the spectral imprint of neutral 
hydrogen. These first galaxies produce copious amounts of 
Lym an alpha radia tion a nd many have now been detected up to 
z=l (ISteidel et a l. 2000; Matsuda e t al. 2004; S aito etal.ll2006L 
l2008l:l Yang e~. 2009: Ouchi et al.b2009: Lehnert et al. 11201 Ol) . 

From a theoretical point of view, the first gal axies are formed 
in atomic cooling hal os with Tyir > 10^ K (Br omm & Yoshidal 
1201 11: iJohnsm] 1201 Ih . Gas in atomic cooling halos may cool 
through Lyman alpha radiation down to 8000 K. Molecular hy- 
drogen is the only efficient coolant in the absence of metals 
that can cool the gas down to a few 100 K. Trace amounts of 
H2 can be formed through gas phase reactions j Saslaw & Zipov 
119671: iLepp & Shulllll984l:TGalli & Pallalll998h . HD is another 
molecule that is a very efficient coolant around 100 K and may 
cool the gas down close to the CMB temperature. The presence 
of deuterium creates a small dipole moment, which allows ro- 
tational transitions to take place with an appreciable Einstein 
A coefficient. Despite the fact that H2 and HD are the effi- 
cient coolants, they are very fragile to Lyman- Werner UV radi- 



ation (Diikstra et aL"'2008"; 'Latif et al.l[20TT^: IShang et al.l[2QTQl: 
IWolcott-Green & Haiman 201 1). 

The first galaxies were preceded by the very first 
stars, forming in minihalos at redshift ~ 20 (I Abel et al 



2002l||Bromm & LarsonI [20041: lYoshida et al.] l2008l: IClark et al 



2011bllal) . These stars emit UV radiation, photoionize gas 
and photodissociate molecules. Consequently, they influence 
the subsequent format ion of structure in t he surrounding 
intergalacti c me dium (I Johnson et al.l l2007l: iJohnsonI l201ll: 
[Maio et al] l201ll) . It is likely that atomic cooling halos are 
polluted by trace a mounts of meta ls through Pop III su- 
pernova explosions (iGreif et al.l l2010l) . The presence of met- 
als or dust significantly influences the thermal evolution 



of ga s and fr agmentat ion be c omes i nevitable ([B romm et al 
2OOII: [Schneider et al.a I2002L l2003l: iBromm & Loeb 20031 



Omukai et al.l l2005l: ISantoro & Shulll [20061: ISmith et all 12008 
Additional metals may originate from the local presence of stel- 
lar winds. The impact of dust and meta ls in the presence of a 
background UV flux has been studied by [Omukai et al.[ ([2008[) . 

The formation of H2 and HD is significantly boosted via 
grain surface reactions as compared to the gas phase route 
in the presence of trace amounts of dust (i.e., > 10"^ solar) 
( Cazaux & Spaans 2004, 2009). Dust grains also attenuate UV 
radiation. The fate of the halo further depends on the con- 
tent of dust, intensity of the UV radiation field, gas-grain col- 
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lisional co upling, and the forma ti on of molecules under thes e 
conditions (iJappsen et al.ll2009allbl : [Schneider & Omukaill201Ql) . 
At high densities, gas-grain collisions have a significant impact 
on the thermal evolution and the stellar initial mass function 
(ISpaans & Silkll200Ql: IPopcke et al.ll201lh . These are therefore 
the key parameters to study. 

Although the number of observed Lyman alpha emitters 
has increased manifold during the past years, the inception of 
Lyman alph a emitters is not ye t fully understood. In our ear- 
lier studies (iLatif et al.ll201 lblla|) , we focused on metal free ha- 
los. However, some o bservations siiggest the presence of dusty 
Lyman alpha emitters (iTamura et al .'2QQ9: Ono et al.ll201Ql) . The 
presence of dust has important implications for Lyman alpha 
emission. It may absorb Lyman alpha photons and re-em i t them 
in t he continuum at far-infr ared frequencies. iNeufeldl (Il99ll) 
an aiman & SpaansI (Il999l) studied the Lyman alpha radiative 
transfer in a dusty medium. They found that Lyman alpha pho- 
tons face less attenuation than nonresonantly scattered photons 
if the medium is inhomogeneous. 

In this paper, we study the impact of dust on the dynam- 
ical evolution of a protogalactic halo, the formation of binary 
disks, and the implications of dust for the emission of Lyman 
alpha photons. We intend to asses how the presence of dust in- 
fluences the formation of H2 and HD molecules for diff'erent 
strengths of the impinging UV radiation field. To complete this 
task, we have performed state-of-the-art high resolution cosmo- 
logical simulation s by extending our previous chemical model 
(iLatif et al .1120 11 ah to include dust physics. The main ingredi- 
ents of our extended chemical model include the formation of 
H2 and HD molecules on dust grains, extinction of UV radiation 
by dust grains and gas-grain collisional coupling. We study the 
evolution for diff'erent dust-to-gas ratios and various strengths of 
the UV radiation field. 

Our paper is organized in the following way. In the 2nd sec- 
tion, we describe the numerical schemes and simulation setup. 
We briefly summarize our chemical model in the 3rd section. We 
present the results obtained in the 4th section. In the 5th section, 
we discuss our conclusions. 

2. Computational Methods 

Our numerical simulations have b een performed usin g the ex- 
tended version of the FLASH code to ubev et al.ll2009l) . FLASH 
is an adaptive mesh refinement (AMR), modular, parallel, grid 
based code. It can run on massively parallel systems and can be 
used for a wide variety of astrophysical problems. It uses the 
message passing interface (MPI) to achieve portability and scal- 
ability on diff'erent systems. It discretizes the computational do- 
main into nested grid cells by making use of the PARAMESH 
library. It has two exchangeable grids, a uniform grid and 
an oct-tree based adaptive grid. We exploit the AMR tech- 
nique to add resolution in the domain of interest. We capital- 
ize the AMR method and add 15 additional levels of refine- 
ment. Consequently, we obtain an eff'ective resolution of 5 pc. 
We resolve the Jeans length by at least 20 cel ls to make sure 
that all essential processes are resolved properly (iFederrath et al.l 
[2QTTbh. This also ensures the fulfillment of Truelove criterion 
(iTruelove et"aDll997l) . We use an unsplit hydro solver with a 
3rd order piece- wise parabolic (PPM) method for hydrodynamic 
calculations. A multigrid Poisson solver is employed for self- 
gravity computations. The particle mesh method is used for the 
dynamical evolution of a dark matter. 

We perform 3 -dimensional cosmological simulations and 
our computational box has a comoving size of 5 Mpc in each 



dimension. Periodic boundary conditions have been used both 
for the hydrodynamics and gravity. Our simulations start with 
cosmological initial c onditions and make use of the COSMIC 
package developed by iBertschinged (1 19951) to produce Gaussian 
random field initial conditions. We start our simulations at red- 
shift 90 as computed by the COSMIC code, employing 3.2 x 10^ 
particles for dark matter gravity calculations. The dark matter 
particles are distributed according to the initial distribution of 
the baryons (see below). We use constraint realization, available 
in the grafic code (part of the COSMIC package), to select a mas- 
sive halo at the center of a box. We start with an eff'ective grid 
resolution of 512^ in the central 1 Mpc region and set the rest of 
the box to a resolution of 128^ grid cells. In this way, the reso- 
lution in the central Mpc is increased by a factor of 4. Our sim- 
ulations are based on a AC DM cosmology with WMAP 5 -years 
parameters Qm = 0.2581, Hq = 72 km s"^ Mpc-\ Qb = 0.0441, 
assuming a scale invariant power spectrum with erg = 0.8. The 
size of the dark matter halo at z=5.4 is -10 kpc and has a mass 
of ~ 2 X 10^^ Mq. Simulation results are partly anal yzed using 
YT, i .e., a visualization toolkit for astrophysical data (iTurk et aP 
l2Qllh . 

3. Chemical Model 

The gas in the protogalactic halos is heated up to their virial tem- 
peratures where cooling due to Lyman alpha radiation becomes 
important. HI can cool the gas down to 8000 K. In the absence 
of metals, H2 is the only molecule that can cool the gas down 
to a few xlOO K. Trace amounts of H2 can be formed through 
the residual fraction of free electrons left over by the recombi- 
nat ion epoch. H2 can also be formed in s hock heated regions 
( Te gmark et a n il997l: lAnninos et al.1ll997h . The gas phase for- 
mation is dominated by 

H e" ^ H" -h y. (1) 

This leads to H2 formation via the reaction 

H-hH-^H2 + e". (2) 

HD is a molecule whose cooling becomes efficient below 200 K. 
The first excited state of HD can be excited around 150 K. It can 
cool the gas down to about Tcmb = 2.75(l-hz). The main route 
for HD formation is 

H2 +D+ ^ HD-hH+. (3) 

In a previous study (iLatif et al.l l2011ah . we devised a 
chemical model for a gas of primordial composition, con- 
sisting of 36 reactions (23 collisi onal and 1 3 radi ative). 
For details, see the appendix of (iLatif et al.l l2011ah . We 
solve the rate equations of the following 12 species: 
H, H+, He, He+, He++, e", H", H2, H+, D, D+, and HD. 
We include non-equilibrium ionization, photoionization, colli- 
sional ionization, radiative recombination and photodissociation 
of the above mentioned species. We also include a comprehen- 
sive model for cooling and heating, i.e., photo-ionization heat- 
ing, collisional excitation cooling, recombination cooling, col- 
lisional ionization cooling, bremsstrahlung cooling and photo- 
dissociation heating. Our model also comprises cooling due to 
H2 and HD molecules as well as a multi-level model for atomic 
hydrogen line cooling (see appendix: table 2 of iLatif et afl 
(201 la)). We also include the self- shielding of H2 and HD. 

We have extended our previous chemical model to include 
dust physics. For the description of the dust abundance, we adopt 
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the dust-to -gas ratio measu r ed in u nits of the solar value, as sug- 
gested by ISchneider et alJ (1201 ll) . Dust models show that su- 
pemovae can efficiently produce dust CTo dini & Ferrara 2001; 
iNozawa et al] l2003l: iBianchi & Schneided l2007l) . Indeed, sig- 
natures of supernova dust have been observed in the spectra 
of a h igh-redshift quasar (iMaiolino et alJ 120041: iGallerani et alJ 
I2010h . 

The presence o f dust efficiently boosts the for mation of H2 
and HD molecules (ICazaux & Spaansll2004 l2009h for metalHc- 
ities as low as Z/Z© = 10"^. In order to include a recipe for 
the formation of H2 and HD molecules on dust surfaces, we 
adopted the H2 and HP formation rates on dust grains from 
ICazaux & SpaansI (l2009l) , see equations 12 and 13. The forma- 
tion efficiency of H2 and HD on dust grain surfaces depends 
on the gas and dust temperatures as well as on the grain size 
distrib ution and composition We u se the grain size distribu- 
tion of IWeingartner & Draind (1200 ll) . In this study, we further 
focus on carbon (PAHs, AC) and silicate dust grains. We have 
grain sizes of 3.5-300 A. We assume that dust abundance scales 
linearly with overall metallicity and that the total dust grain 
cross section is equal to the the Milky Way for solar metallicity 
(ICazaux & SpaansI [20091) . In this way, we take into account the 
contribution of smaller grain sizes which may further enhance 
the H2 and HD formation rates due to an increase in surface area. 
The dust tem perature in th e presence of a background UV flux is 
computed as (iHoUenbacFet al 199TI:lMeiierink & Spaansl2005h 

Td = {8.9 X 10"^ VoGolO"^-^^^ T^j^B 3.4 x 10"^ (4) 
[0.42 - log(3.5 X lO-^iooTo)] x noo X/^f ^ 

where vq = 2.65 x 10^^ s"\ Go is the far-ul traviolet (FUV) 
flux in units of the equivalent Habing flux (iHabingl [l968i) . 
To = 12.2Gq^ K is the equilibrium dust temperature and tioo = 
2.7 X IOVT^GoZ/Zo is the optical depth at 100 yum. The dust 
temperature is computed assuming that each dust grain is illu- 
minated by attenuated FUV flux, cosmic microwave background 
radiation ( CMB), and the infrared radiation field from warm dust 
emission (iHollenbach et al.l [19911) . If gas and dust grains have 
diff'erent temperatures, they can exchange heat through colli- 
sions. We compute the gas-grain col lisional heating and cooling 
rates as given in the equation below (iMeijerink & Spaansll2005l) 



r.n=1.2xlO-„^(Z/Ze)(^) (^) (5) 
x[1.0-0.8exp(-75/Tk)](Td-Tk) 



ate molecules. We have modeled th e background UV flux in the 
following way (lOmukai et al.ll2008h 

J^" = J2ilO-2HB,(T,)/B,H]f(v) erg cm-^sr-^s'^Hz-^ (6) 

where J21 is the background UV radiation field intensity below 
the Lyman limit and T=, is the color temperature of a star. We 
assume a black body radiation spectrum with T., = 10"^ K. f(y) 
is one for non-ionizing radiation and is equal to fesc for ioniz- 
ing radiation. ByH is the black body spectrum at the Lyman limit 
(J21 = 10"^^ erg/cm^/ s/Hz/sr). By is the stellar radiation spec- 
trum (i.e., black body spectrum) as given by 

By (T) = ^ . (7) 

c^ (exp(hy/kT) - 1) 

Here h is Planck's constant, v is the frequency of radiation, c is 
the speed of light, k is Boltzmann's constant and T is the radia- 
tion temperature. Jih = J21 x fesc is the ionizing flux, where fesc 
is again the escape fraction of ionizing radiation. These ionizing 
photons are attenuated while propagating through intergalactic 
medium. We have taken this efl'ect into account in an approxi- 
mate manner by computing the optical depth for each frequency 
bin in the following manner 

Ti(yj) = i:iO-i(yj)xniX.lj, (8) 

where i stands for the specie and j denotes the frequency bin, 
(Ti is the photoionization cross section of the specie, n is the 
gas number density and Aj is the local Jeans length. The radi- 
ation flux also becomes attenuated by dust grains. These grains 
both absorb and scatter the radiation. The dust opacity strongly 
depends on the fre quency of the UV flux . We use a fit to the 
extinction curve by iRoberge et al.l (Il99ll) to compute the dust 
opacity as given by 

Tdust(vi) = nH(i)o-d(i)Ay(i)^j x Z/Z© (9) 

here nu is the number density of hydrogen, cr^ is the eff'ective 
dust cross section, and Ay is the dust extinction coefficient. Total 
attenuation of the radiation field is given by 

r\y) = r-\y) X exp(-T(yO) X exp(-Tdust(n)). (10) 

The stellar sources formed inside the halos will produce a lo- 
cal UV flux. We have also modeled this flux and studied its influ- 
ence on the collapse and chemistry of the gas physics. As these 
sources are formed in the high density regions, we enhance the 
background UV flux by factors of 100 and 1000 above densities 
of 1 cm"^. 



Td and Tk are dust and gas temperatures, respectively, n is the 
gas number density, and amin is the mini mum dust grain size. 
We s et the minimum grain size to 3.5 A (IMeijerink & SpaansI 
[2OQ5I) . Grain surface reactions are often more efl'ective than their 
gas-phase counterparts. We have included H2, HD, H, D and He 
formation on the surfaces of dust grains in o ur chemical model. 
The ra tes f or these reactions are tak en from iGlover & JappsenI 
(l2007h and ICazaux & SpaansI (l2009l) . and are Hsted in table [B 
We have also included the H2 collisional dissociation cooling, 
photoelectric heating by UV irradiated dust grains, gas-phase H2 
formation heating, and heati ng by H2 formation on d ust grains. 
These are listed in table 7 of iGlover & JappsenI (l2007l) . 

During the epoch of reionization, a background UV flux 
is produced by stellar populations. The background UV flux 
can photoheat and photoionize the gas, and can photodissoci- 



4. Results 

4.1. Global dynamics 

We have performed a grid of simulations for different strengths 
of external and internal UV ffuxes as well as for two different 
dust-to-gas ratios as shown in table [21 We start our simulations 
at redshift 90 with cosmological initial conditions. Density per- 
turbations decouple from the Hubble flow and begin to collapse 
through gravitational instability at redshift 20. The gas falls into 
the dark matter potentials and gets shock heated. The virializa- 
tion processes partly transform the gravitational binding energy 
into kinetic energy of the gas and dark matter. Part of the gravi- 
tational potential energy goes into thermal energy of the gas and 
is radiated away. We also see the merging of smaller clumps into 



3 



Table 1. 



Latif et aL: The chemistry and dynamics during the formation of high-redshift protogalaxies 
CoUisional and radiative rates 



No. 


Reaction 


Rate 


Reference 


1 


H2 formation rate on dust grains 


see reference 


CS09 


2 


HD formation rate on dust grains 


see reference 


CS09 


3 


H+ + e-^ H 


see reference 


GJ07 


4 


D+ + e"^ D 


see reference 


GJ07 


5 


He+ + e-^ He 


see reference 


GJ07 



Glover and Jappsen (2007);CS09 :Cazaux and Spaans 2009. 




Fig. 1. Time sequence of density snapshots at different redshifts. Redshift and distance scales are mentioned in each panel, and 
shown for model A. The values of density shown in the colorbar are in comoving units [g/cm^]. The merging of clumps is shown 
inside the black circle. 




19.0 19.5 20.0 20.5 21.0 21.5 22.0 19.0 19.5 20.0 20.5 21.0 21.5 22.0 

Log 1 Radius [cm] Log 1 Radius [cm] 

Fig. 2. Average density radial profiles for diff'erent strengths of the background UV field and dust-to-gas ratios at z=5.4. The values 
of the background UV flux and metallicities are shown inside each figure. The right hand side panel is for an internal UV radiation 
flux (model E, F, G, H) in addition to a background UV flux of J21 = 1. The value of the escape fraction for ionizing radiation is 
10% for both panels. 



bigger ones. The occurrence of a major merger is shown in fig- 
ure [T] This process of merging helps in the transfer of angular 
momentum, which is necessary for th e gas to settle in the grav- 
itational potential (iMayer et al .1120101) . We switch on the radia- 
tive backgrounds at z = 8, where we expect reionization to be 



already in its final phase. The temperature of the gas is increased 
above 10^ K due to photoheating, depending on the strength 
of the radiation field. At these temperatures, the recombination 
cooling kicks in and the gas begins to recombine. The gas in 
the filaments is relatively denser, becomes self- shielded and falls 
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Fig. 3. Average radial profiles for diff'erent background UV fields and dust-to-gas ratios at z=5A. The upper left panel of this figure 
shows the temperature radial profile of the halo. The HII abundance radial profile for the halo is depicted in the upper right panel. 
H2 abundance is shown in the lower left panel. The lower right panel shows the HD radial profile of the halo. The escape fraction 
of ionizing radiation is 10% for all the panels shown in the figure. 

Table 2. Estimates of dust continuum emission for diff'erent UV radiation fields 



Model 


Background UV flux 


Internal flux 


Dust-to-gas ratio 


Dust temperature 




h 1 [erg / cm^ / s /Hz / sr] 


hi [erg/cm^ /s /Hz/ sr] 


[Z/Zo] 


[K] 


A 


1 





1.0 X 10-' 


18.4 


B 


1 





1.0 X 10-4 


18.63 


C 


10 





1.0 X 10-2 


18.6 


D 


10 





1.0 X 10-4 


18.7 


E 


1 


100 


1.0 X 10-2 


22.4 


F 


1 


100 


1.0 X 10-4 


22.5 


G 


1 


1000 


1.0 X 10-2 


33.2 


H 


1 


1000 


1.0 X 10-4 


33.6 



into the center of the galaxy th rough cold streams (iHaiman et al.l 
I2OOOI: iDiikstra & Loebll2009h . These cold streams have typical 
temperatures of 10"^ K and number densities on the order of 
0.01 - 1 cm"^. The average density radial profiles for diff'erent 
strengths of the background UV fields are shown in figure[2j The 
variations in the density radial profiles are due to the substruc- 
ture in the halo and vary with the strength of the radiation field 
and dust-to-gas ratio. Before discussing the disk morphologies 
and dynamics in more detail, we will explore the impact of dust 
and radiation backgrounds in section 4.2. 



4.2. Chemical and thermal evolution 

The average temperature radial profiles for two diff'erent dust- 
to-gas ratios and various strengths of a background UV radia- 
tion field are shown in the upper left panel of figure [S] It can be 
seen that gas is photoheated by a UV radiation flux in the en- 
velope of a halo and consequently its temperature is raised to 
above 10^ K. It can also be seen that the maximum temperature 
is 10^ K for ionizing flux (J21 xfesc) of 1 in units of J21. The cool- 
ing due to helium lines becomes important at > 5 x 10"^ K and 
thus keeps the gas temperature at 10^ K. Recombination cooling 
becomes eflTective above 10"^"^ g/cm^ and cools the gas down 
to 10"^ K. It is found that for model A, gas cools eflSciently due 
to the enhanced fraction of H2 and HD molecules in the pres- 



5 



Latif et aL: The chemistry and dynamics during the formation of high-redshift protogalaxies 



4' 



I 



X-axis [2 kpc] 



I 



Fig. 4. Density weighted projection of gas density for different background UV fields and dust-to-gas ratios at redshift 5.4. The top 
left panel shows the density for model A. The density structure for model B is shown in the upper right panel. The bottom left panel 
shows the density structure for model C while the bottom right panel shows the density structure for model D. The escape fraction 
of ionizing radiation is 10% for all the panels. The values of density shown in the colorbar are in comoving units [g/cm^]. 



ence of dust grains. The H2 and HD cooling begins to become 
effective around 10"^^ g/cm^, where gas is cooled down to -150 
K. For model B, the temperature in the center of the halo is in- 
creased to approximately 250 K. Similarly, for models C and D 
(i.e., J21 = 10) the temperature is higher in the center of a halo 
due to the photodissociation of H2 and HD molecules. 

The temperature evolution is significantly different for 
higher and lower dust-to-gas ratios. The ionization degree of the 
gas is shown in the top right panel of figure [3l The ionization 
fraction of HII is higher in the surroundings of the halo, de- 
pending on the intensity of ionizing radiation, and goes down 
towards the center of the halo. The overall degree of ionization is 
higher for a stronger radiation field and lower for a weaker radi- 
ation field. At radii < 10^^-^ cm, the ionization fraction decreases 
sharply as gas becomes neutral and is efficiently shielded from 
the ionizing rad iation flux. The re sults are in agreement with our 
previous study (iLatif et al. 11201 lab . 

The bottom left panel of figure [3] depicts the H2 radial pro- 
files. It is found that the molecular hydrogen abundance is very 
low in the envelope of the halo as gas is highly ionized. There 
is a sharp rise in H2 abundance at radii between 10^^-^ cm 
and 10^^ cm. This is due to the substructure in the halo where 
gas self- shields against the background UV radiation field and 
molecular hydrogen survives. The H2 abundance sharply in- 



creases towards the center of the halo as its photodissociation 
becomes ineffective due to further self- shielding. For a weaker 
radiation field and higher dust-to-gas ratio (i.e., model A), the 
H2 fraction is higher compared to the stronger radiation field 
and lower dust-to-gas ratio cases. The H2 abundance is signifi- 
cantly boosted in the core of the halo due to the formation of H2 
molecules on dust grains at high gas densities. This increase in 
H2 abundance in the presence of dust grains is consistent with 
earlier studies of Cazaux & Spaans ( 2009, 2004). Similar be- 
havior is found for HD molecules, as shown in the bottom right 
panel of figure [S] The HD abundance is lower for a stronger UV 
background flux and higher for a weaker flux. A maximum HD 
abundance of 10""^-^ is found for model A. The HD abundance 
remains lower than the critical value (i.e., 10"^) for J21 > 0.1 
in agreement with prev ious studies (IWolcott-Green & HaimanI 



m agreement witn prey i' 
l20ritlL^tifera^ 



4.3. Disk formation and morpiiology 

The properties of simulations for different strengths of back- 
ground UV flux and dust-to-gas ratios at redshift 5.4 are illus- 
trated in figure |4l It is found that binary disks are formed for 
model A, as shown in the top left panel of figure IH in differ- 
ent subhalos. We note that the upper disk is formed first and 
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Fig. 5. Average radial profiles of disks for diff'erent background UV fields and dust-to-gas ratios at z=5.4. The left panels of this 
figure show radial profiles for disk 1 and the right panels show the radial profiles for disk 2. The top two panels show the ratio 
of rotational to orbital velocity for the disks. The middle panels show the ratio of turbulent to thermal energy for each disk. The 
Toomre parameter Q for each disk is depicted in the bottom panel. The escape fraction of ionizing radiation is 10% for all panels. 



the second disk appears a few megayears after the formation of 
first. Each disk has spiral arms and gas masses of 4 x 10^ M© 
and 2 X 10^ M©, respectively, and sizes of ~ 1 kpc. We see 
indications for the formation of an additional disk during later 
stages of the evolution. For the same radiation field, but with 
lower dust-to-gas ratio (i.e., model B), the halo fragments again 
and forms binary disks. The morphology of the halo is shown 
in the top right panel of figure IH and is very diff'erent from the 
previous case. The temperature of the halo is higher due to the 
lower dust-to-gas ratio. The gas masses of the lower and up- 
per disks are 2.3 x 10^ M© and 1.4 x 10^ M©, respectively. We 
found that the lower disk is formed first (opposite to the pre- 



vious case) and the second disk is formed after 6 megayears. 
The bottom left panel of figure IH shows the structure of a halo 
for model C. It can be seen that again binary disks are formed. 
The formation process is diff'erent, the top disk is formed first 
while the second disk is formed through the merging of smaller 
clumps. Moreover, the formation of the bottom disk is delayed 
compared to the lower radiation background UV field case and 
is formed 19 megayears after the formation of the first disk. The 
overall collapse of the halo is also delayed due to the stronger 
radiation field. The gas masses of the lower and upper disks are 
2.0x 10^ M© and 3x 10^ M©, respectively. The bottom right panel 
of figurelH shows the state of the halo for model D. Similar to the 
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Fig. 6. Average radial profiles for diff'erent internal UV field strengths and dust-to-gas ratios at z=5.4. The upper left panel of this 
figure shows the temperature radial profile of the halo. The HII abundance radial profile for the halo is depicted in the upper right 
panel. The lower left panel shows the H2 abundance. The HD radial profile of the halo is depicted in the lower right panel. The 
escape fraction of ionizing radiation is 10%. 



other cases, two disks are formed. The upper disk is formed first 
and the second disk is formed after 20 megayears. This delay 
is due to the enhanced strength of the radiation field and lower 
dust-to-gas ratio. Gas masses of the lower and upper disks are 
2.2 X 10 Mo, 1.7 X 10^ Mo, respectively. The binary disk forma- 
tion in our simulations is due to the presence of subhalos after a 
major merger. Employing diff'erent initial conditions may result 
in single or multiple disks. 



4.4. Disk dynamics and Stability 

We have computed the circular velocity (Vc = ^/GM/R, M is 
the disk mass) and rotational velocity (Vrot = L/R, where L is 
the specific angular momentum of the disk) for each disk. The 
ratio of rotational to circular velocity for disk 1 (the massive 
disk in each case) is shown in the top left panel of figure [51 
while for disk 2 (the less massive one in each case) is shown 
in the top right panel of figure [S] The rotational velocity remains 
lower than the circular velocity which means that the collapse 
is never halted by rotational support. In the very central region, 
there is a central massive clump with no spiral arm structure, and 
with very little rotational support. This is because the turbulent 
velocities are much larger than the rotational velocity, leading 
to a spherical rather than a disk- type structure. At larger scales 
> 10^^-^ cm , the rotational velocity is increased and we see a 
disk with spiral arms. We found that both disks are co-rotating, 
independent of metallicity. The sign change in the rotational ve- 



locity is a result of the two spiral arms, as the arms between 
the disks then move in opposite directions. We found that disk 2 
in each case has higher rotational support than disk 1 . We also 
computed the turbulent energy for the system and compared it 
with the thermal energy. The ratio of turbulent to thermal en- 
ergy is plotted for both disks in the middle panel of figure [3 It 
is found that turbulence is highly supersonic and turbulent pres- 
sure is more important than the thermal pressure. This is due 
to the high virial temperature, which is refl ected in the kinetic 
energy of the gas. Our results are similar to I Aykutalp & SpaansI 
(2011). To investigate the stability of each disk, we computed 
the Toomre instability parameter Q, given by 



where Q is the Keplerian orbital frequency, Eg is the surface den- 
sity of the disk, Cg = -h v^^^, Cs is the isothermal sound speed 

and Vrms is the root mean square velocity Vrms = -^^^T^^^^tang and 
G is the gravitational constant. The bottom left panel of figure [5] 
shows the spatial evolution of the Q parameter for disk 1 . It is 
found that only for R < 10^^-^ cm the value of Q is larger than 
one and the disk is stable. The sharp increase towards the center 
of the disks is due to the increase in epicyclic frequency. The 
Toomre instability parameter Q for disk 2 is shown in the bot- 
tom right panel of figure [51 The value of the Toomre parameter 
is not well defined on large scales (i.e., > 10^^ cm), as there is 
no well defined disk and the overall spatial evolution of mass is 
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Fig. 7. Density weighted projections of gas density for various internal UV field strengths and dust-to-gas ratios at redshift 5.4. The 
top left panel show the density for model E. The density structure for model F is depicted in the upper right panel. The left panel 
shows the density structure for model G. The state of a halo for model H is shown in the right panel. The escape fraction of ionizing 
radiation is 10% for all panels. The values of density shown in the colorbar are in comoving units [g/cm^]. 





Fig. 8. Density weighted projections of temperature for diff'erent internal UV field strengths and dust-to-gas ratios at redshift 5.4. 
The left panel shows the temperature for model G. The temperature for model H is depicted in the right panel. The panels show 
temperatures corresponding to the density panels of figure|71 This shows that thermal instability leads to the formation of gas clumps 
at these temperatures. The values of temperature are shown in the colorbar. 



much larger than the Jeans mass. We also computed the specific 
angular momentum spatial evolution. We found that angular mo- 
mentum increases with radius according to expectation. We have 
a specific angular momentum on the order of 10^^ cm^/s and our 
results are again in agreement with Aykutalp & Spaans (2011). 
Our stability analysis shows that disks are stable and are sup- 
ported by turbulent pressure. In spite of the global stability, the 



high turbulent Mach numbers considered will locally produce 
strong overdensities in shocks, and may lead to a star forma- 



tion mode regulated by turbulence ([Mac Low & Klessenll2004l: 
iKrumholz & T^l2007fc iBraun & Schmidtll201lh . The final fate 



of the disks depends on the ambient conditions. They may merge 
with each other and lead to the formation of a protogalaxy. They 
may form stars if they become Toomre unstable. 
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Fig. 9. Average radial profiles of the disks for diff'erent strengths of internal UV fields and dust-to-gas ratios at z=5.4. The left panel 
of this figure shows radial profiles for disk 1 and right panels shows the radial profiles for disk 2. The top two panels show the ratio 
of rotational to orbital velocity for the disks. The Toomre parameter Q for each disk is depicted in the bottom panels. The value of 
an escape fraction for ionizing radiation is 10% for both panels. 



4.5. The implications of internal radiation fields 

We also performed simulations with a given strength (i.e., J21 = 
1) of the background UV radiation field, but for an enhanced in- 
ternal radiation flux of 100 and 1000 times the background at 
densities > 1 cm"^ due to the formation of stellar sources inside 
the halo. Figure [6] shows the radial profiles for temperature, ion- 
ization degree of the gas, abundances of H2 and HD molecules. 
It is found that the temperature of the gas is relatively higher as 
compared to the cases with no internal radiation sources. The 
degree of ionization is also enhanced. The gas still cools and 
recombines inside the HIT regions, which is why the degree of 
ionization declines towards the center of a halo. The H2 and HD 
abundances are shown in the bottom panel of figure[6l It is found 
that due to the enhanced strength of the internal UV flux, pho- 
todissociation of molecules becomes efl'ective and H2 and HD 
are almost two orders of magnitude lower than in the no internal 
radiation case. The state of the simulations for internal radiation 
sources at redshift 5.4 is shown in the figure |7] It is found that 
the collapse of a halo is delayed for these cases. We find that for 
model E, gas is still able to cool down to a few hundred Kelvin, 
while for the lower dust-to-gas ratio case (i.e., model F) the tem- 
perature increases up to 1000 K. A binary system is again formed 
irrespective of the dust-to-gas ratio. The gas masses of disk 1 
and 2 are 1.9 x 10^ Mo and 1.7 x 10^ M© for Z/Zq = 10"^ 
and 1.6 x 10^ M© and 8 x 10^ for Z/Z© = 10"^ Despite the 
higher temperatures in the core of a halo due to photodissoci- 
ation of molecules, surprisingly, fragmentation is enhanced for 



J21 = 1000 irrespective of the dust-to-gas ratio. This is due to 
cooling instabilities because the temperature lies in the unstable 
regime of 2000 to 5000 K. The formation and dissociation of H2 
is very sensitive to the above-mentioned temperature range. As 
the temperature increases above 3700 K, the H2 survival time 
sharply decreases and the cooling time scale increases. On the 
other hand, if the temperature is slightly decreased, this leads to 
longer H2 survival times (see figure 4 of Oh&Haiman (20(3)). 
This process initiates thermal instabilities. The gas masses of the 
fragments are on the order of a few 10^ M©. Figure [8] shows the 
temperature corresponding to the density projections in figure [Tj 
(i.e., models G and H). It can be seen that in the high density re- 
gions fine-grained thermal substructure occurs, which indicates 
small pockets of gas that quickly cool down to low temperatures. 
This supports our hypothesis of thermal instabilities leading to 
the local cooling and collapse of gas. 

We also computed the ratio of rotational to circular veloc- 
ity of the disks for internal radiation field cases (i.e., models 
E, F, G and H). These are shown in the top panels of figure 
The ratio remains less than unity as in previous cases. We 
find that for a lower dust-to-gas ratio, the disks are now counter- 
rotating while they are co-rotating for higher dust-to-gas ratio 
case. Dyn amical evolution is faster in high dust-to-gas ratio envi- 
ronments (lAykutalp & Spaansll2Qlll) . Therefore, the lower-mass 
disk falls behind in its evolution compared to the higher-mass 
one for lower dust-to-gas ratio. Subsequently, the lower-mass 
disk is forced to accrete gas from the trailing spiral arm of the 
higher-mass disk and counter-rotation ensues. We also computed 
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the turbulent and thermal energy of the system and found that 
turbulence is again highly supersonic and the disks are supported 
by turbulent pressure. The Toomre stability parameter Q is larger 
than unity in the central region (i.e., < 10^^^ cm), where we see 
a well-defined disk. The disk is therefore expected to be stable, 
similar to previous cases, as shown in the bottom panels of figure 

il 

4.6. Lyman Alpha emission from high-redshift protogalaxies 

We estimate the Lyman alpha flux emanating from the galaxy 
forming halo. Our results for the Lyman alpha emission in the 
presence of dust are shown in figure [TOl We found that above 
columns of 10^^ cm^ and optical depths of r > 10^ trapping of 
Lyman alpha photons becomes eff' ective, co nsistent with previ- 
ous studies ( Spaans & Silk 2006; Latif et al.t 2011a.b.c). That is 
why Lyman alpha emission sharply declines towards the center 
of the halo. It can be seen that the Lyman alpha flux emerges 
from the envelope of the halo. The value of the flux depends 
upon the strength of the external UV radiation field and the dust- 
to-gas ratio. The local variations in the flux are due to the dumpi- 
ness of the medium. We assume a constant dust-to-gas ratio, 
leading to a highly inhomogeneous, clumpy interstellar medium 
as a resu lt of the high Mach numbers in our simulation (see 
models of lNeufeldll99llHaiman & Spaansll99^IShimizu et alJ 
1201 Ih . As dust sits in the center of the dense clumps, Lyman al- 
pha photons scatter ofl' these clumps and are able to escape the 
halo due to large scattering cross sections. This depends on the 
clump covering factor, i.e., average number of gas clouds along 
the line of sight. We compute the column density distribution 
and use it to determine the clump covering factor. We use the 
model o fE aiman & SpaansI (1 19991) to compute the HI L yman al- 
pha escape fraction. Our findings are in agreement with lNeufeldl 
0991); Haiman & Spaans (1999); Shimizu et al. (201l|). 

The presence of a background UV flux excites the emission 
of Lyman alpha photons and enhances its luminosity. Moreover, 
flux predominantly arises from collisional excitation of atomic 
hydrogen and the contribution from recombination line emission 
is a few orders of magnitude lower. For further details on the ori- 
gin and sou rce of Lyman alpha emission see our previous studies 
(lLatifetal.ll20TTbllah . Our results also agree with previous stud- 
ies that emission of Lyman alpha is extended and emerges from 
the warm gas in the surro i indings of the galaxy (Diikstra et alj 
2006allbl: iDiikstra & Loebl l2009l: iFaucher-Giguere et al. 201( ^ 



Latif et al.ll2011ah . We found that a flux on the order of 10"^^ _ 



10"^^ erg/cm^/s emanates from the the protogalactic halo at red- 
shift 5.4. The redshifted wavelength of Lyman alpha photons 
at redshift 5.4 is 0.78 micron. This should be detectable with 
present-day telescopes like Subaru Q and VLTB Further details 
on the issue of line trapping of Lyman alpha photons and the role 
of backgroun d UV flux in the e mission of Lyman alpha photons 
are given by ISchleicher et al.l (l2010bb and Latif et al. ( 2011a). 
Lyman alpha emission, apart from dust content, also depends on 
the velocity structure of a halo. The variation in Doppler velocity 
shift does no t strongly influence the es cape fraction of Lyman al- 
pha photons (iHaiman & SpaansI 19991) . We also checked the con- 
sistency of intrinsic Lyman alpha luminosity with the adopted 
internal UV flux. We computed flux from intrinsic Lyman alpha 
luminosity by normalizing it with Lyman alpha line width and 
found it consistent with the adopted UV flux. In the presence of 
dust, the contribution of the internal stellar feedback, on top of 



the background UV field, can vary strongly. For Z/Z© = 10"^, 
UV continuum optical depths range from 10"^ to almost unity, 
and are a factor 100 times larger for Z/Z© = 10"^. In our simula- 
tions, column densities ra nge from 10 ^ ^ to 10 ^"^ cm"^, a similarly 
wide range as found by IWada et aP (l2009l) . For the dust- free 
case, we expect internal stellar feedback to be quite efl'ective, 
and to partly photo-dissociate H2. However, in the presence of a 
strong background UV field, the role of the intrinsic luminosity 
again diminishes. 



4. 7. Dust emission from high-redshift protogalaxies 

We further compute the dust temperatures for diff'erent UV back- 
ground fluxes and dust-to-gas ratios for the cases listed in table 
[21 We found that for lower background UV fluxes (i.e., J21 = 1), 
the dust temperature remains close to the Tcmb while for higher 
fluxes (i.e., J21 = 100) it is well above Tcmb- We have a max- 
imum dust temperature of 34 K for the model H in our case. 
Our results are in agre ement with observed dust temperatur es of 
high-redshift quasars (lOmont et al.ll200TI: lYang et al.ll201 Ih . We 
found that for the given dust-to-gas ratios and the halo masses 
considered here (dust masses, 10^ M© for Z/Z© = 10"^ and 
10^ M© for Z/Z© = 10"^), the dust continuum emission is far be- 
low the detection limits of the Atacama Large Millimeter Array 
(ALMAfl at z=5.4. However, ALMA will be able to explore the 
dust content in larger systems or in the presence of higher dust- 
to-gas ratios. Our calculations show that a 3 sigma detection is 
possible with ALMA band 9 in five ho urs of integrat i on if d ust 
masses (for calculation see equation 1 of lOmont et al.l (l200lh ) of 
6.75 X 10^ M©, 3.27 x 10^ M©, and 1.54 x 10^ M© are present 
(assuming dust temperatures of 20 K, 30 K, and 50 K, respec- 
tively). One requires roughly solar metallicity for a ~ 10^^ M© 
protogalactic halo to be detectable at z > 5. 



5. Discussion and Conclusions 

In all, we have performed 9 cosmological simulations to study 
the dynamics of disk formation in protogalactic halos under the 
influence of difl'erent UV radiation field strengths and dust-to- 
gas ratios. We used the adaptive mesh refinement code FLASH 
to carry out this task and coupled it with a chemical network. Our 
chemical model includes photoionization, collisional excitation, 
collisional ionization, radiative recombination, photodissocia- 
tion processes and the formation of molecules on dust grains as 
well as a multi-level treatment of the hydrogen atom. We found 
that gas in the atomic cooling halo is photoionized by a back- 
ground UV flux and that photoheating raises the gas temperature 
to above 10^ K. At densities > 10"^"^ g/cm"^, gas cools, recom- 
bines and falls into the center of a halo through cold streams of 
gas. The later thermal evolution of a halo depends on the strength 
of a background UV flux and the dust-to-gas ratio. The forma- 
tion of H2 and HD molecules is significantly enhanced (by two 
orders of magnitude) in the presence of dust grains as compared 
to gas phase reactions. We found that the formation of molecules 
is more eff'ective for higher dust- to gas ratio (i.e., Z/Z© = 10"^) 
as compared to lower dust-to-gas ratio (i.e., Z/Z© = 10"^). The 
minimum temperature in the core of a halo is ~ 200 K for model 
A. This value is higher (800 K) for lower dust-to-gas ratio and 
stronger UV radiation flux. We conclude that the thermal evolu- 
tion of a protogalactic halo is strongly influenced by the presence 
of dust and a UV radiation field. 



^ http://www.naoj.org/Observing/Instruments/IRCS/camera/filters.html 3 http://www.almaobservatory.org/en/science-with-alma/evolved- 
^ http://www.eso.org/public/teles-instr/vlt.html stars- Stardust 
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Fig. 10. Lyman alpha emission for different magnitudes of the background UV field and dust-to-gas ratios. The right hand side panel 
shows the flux for difl'erent background UV radiation fields. The left panel shows the Lyman alpha flux for an internal UV radiation 
field in addition to a background UV flux J21 = L The value of the escape fraction for ionizing radiation is 10% for both panels. 



Our results show that binary disks are formed in the proto- 
galactic halo at redshift 5.4 irrespective of the background UV 
radiation field strength and dust-to-gas ratio in our simulations. 
However, the morphology and the formation time of disks does 
strongly depend on the dust-to-gas ratio and strength of a back- 
ground UV flux. We found that disks are stable as the value 
of Toomre instability parameter Q is greater than one at radii 
< 10^^-^ cm. We also found that the turbulence in the halo is 
highly supersonic due to the higher virial temperature of the 
halo. Moreover, disks are more supported by turbulent pressure 
than thermal pressure. It may happen that disks merge with one 
another at lower redshifts. There is also a possibility that more 
than 2 disks are formed in some cases, which will lead to the 
formation of a multiple disk system. We could not follow the 
evolution of a halo to a lower redshift than 5.4 as this becomes 
a computationally too expensive task. Nevertheless, we have ex- 
plored an important phase of disk formation, which may be ob- 
served with future planned surveys. 

We also performed a set of simulations with difl'erent 
Gaussian random field initial conditions to determine whether 
the origin of binary disk formation comes from cosmological 
initial conditions. We found that a single disk is formed in the 
center of a halo instead of binary disks for a diff'erent random 
seed. We thus think that the origin of binary disk formation is 
associated with the underlying dynamics of the dark matter in 
the halo. In an upcoming paper, we will address the issue of how 
often these binary disks form in atomic cooling halos. 

We performed additional simulations for the cases where the 
background UV flux is enhanced by factors of 100 and 1000 
due to the formation of radiation sources like stars inside the 
halo. It was found that the collapse of a halo is delayed in HII 
regions formed by the internal radiation sources and the pho- 
todissociation of molecules (H2 and HD) is strongly enhanced 
in this case. We found that despite the higher temperatures in 
the halo for an internal radiation flux of J21 = 100, in addi- 
tion to a background UV flux of J21 = 1, again binary disks 
are formed for both Z/Z© = 10"^ and Z/Z© = 10""^. For a 
higher strength of the internal radiation flux, J21 = 1000, the 
temperature of the halo in the center is increased to a few thou- 
sand Kelvin. We found that for this case cooling instabilities 
lead to the formation of multiple clumps for both dust-to-gas 
ratios, Z/Z© = 10"^ and Z/Z© = lO""^ at redshift 5.4. Due 
to significant local overdensities in the high-Mach turbulence, 
we expect a star formation mode that is regulated by the turbu- 



lent densitv PDF (iMac Low & KlessenI 120041: [Krumholz & TanI 
I2OO7I: iBraun & Schmidtll2012h . The further fate of these clumps 
depends on how they transfer angular momentum. They may 
merge with each other and form a disk if they lose angular mo- 
mentum. It may also happen that they keep cooling, and collaps- 
ing, and fragment further. 

We estimated the Lyman alpha emission from protogalac- 
tic halos for diff'erent strengths of UV radiation fields and dust- 
to-gas ratios. We found that even in the presence of dust, for 
Z/Z© = 10"^ and Z/Z© = 10"^, Lyman alpha photons scatter off' 
dusty gas clouds and an inhomogeneous dust distribution leads 
to the efficient escape of Lyman a lpha photons. Our r esults are 
in agreement with the models of iHaiman & SpaansI (Il999h . It 
was found that a ffux of the order of 10"^^ - 10"^^ erg/cm^/ s 
emerges from the envelope of the halo due to effective line trap- 
ping of Lyman alpha ph otons in the center of a halo, consis- 
tent wi t h previous studies JSpaans & SiiEl2006l: ISchleicher et al.l 
l2QlQbl: iLatif et al.ll201 lall^T Such a ffux should be observable 
with present-day telescopes like Subaru and the VLT. Our results 
provide a conservative lower limit to the Lyman alpha ffux as we 
have computed the ffux for minimum velocity gradient along the 
line of sight. Radiative transfer of Lyman alpha photons strongly 
depends on the amount of structure in the halo as well as on 
velocity gradients. Therefore, cosmological simulations coupled 
with Lyman alpha radiative transfer should be performed in the 
future to obtain more robust results. 

Finally, we computed the dust emission from the protogalac- 
tic halo in our simulations. We found that due to the low dust 
masses such systems cannot be observed with ALMA at z=5.4. 
Only galaxies with higher metal content and high dust tempera- 
tures can be probed. However, we also found that dust masses of 
few xlO^ M© yield ~ 0.1 mJy of ffux and thus a 3-sigma detec- 
tion in 5 hours. Given that our protogalactic halos are likely to 
build up a solar content of metals within 1-2 Gyr, they should be 
detectable for z < 5. 

The formation of disks may also have signiffcant implica- 
tions for magnetic fields in the first galaxies. In the absence 
of rotation, the small-scale dynamo is expected to produce tan- 
gled magnetic fields, cl o se to equipartition with t urbulent energy 
dArshakian et al. 1 12QQ9I: ISchleicher etaP l201Qah . This process 
was shown to be efficient for a large range of turbulence models 
and Mach numbers (Federrath et al. 201 la; Schober et al. 2 0T2I) . 
The presence of rotating disks, on the other hand, gives rise to 
large-scale dynamo effects amplifying the mean magnetic field 
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(iBrandenburg & Subramanianll2005h . The presence of magnetic 
fields in high-redshift galaxies will be probed in the future with 
the SKA0 
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